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Three-Dimensional Evolution of the Parker Instability under a Uniform Gravity 

Jongsoo S. S. Hong^, Dongsu Ryu^, and T. W. Jones'^ 

ABSTRACT 

Using an isothermal MHD code, we have performed three-dimensional, high- 
resolution simulations of the Parker instability. The initial equilibrium system is 
composed of exponentially-decreasing isothermal gas and magnetic field (along the 
azimuthal direction) under a uniform gravity. The evolution of the instability can be 
divided into three phases: linear, nonlinear, and relaxed. During the linear phase, the 
perturbations grow exponentially with a preferred scale along the azimuthal direction 
but with smallest possible scale along the radial direction, as predicted from linear 
analyses. During the nonlinear phase, the growth of the instability is saturated and 
flow motion becomes chaotic. Magnetic reconnection occurs, which allows gas to cross 
field lines. This, in turn, results in the redistribution of gas and magnetic field. The 
system approaches a new equilibrium in the relaxed phase, which is different from 
the one seen in two-dimensional works. The structures formed during the evolution 
are sheet-like or filamentary, whose shortest dimension is radial. Their maximum 
density enhancement factor relative to the initial value is less than 2. Since the radial 
dimension is too small and the density enhancement is too low, it is difficult to regard 
the Parker instability alone as a viable mechanism for the formation of giant molecular 
clouds. 


Subject headings: instabilities — ISM: clouds — ISM: magnetic fields — ISM: structure 
— MHD 
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1. INTRODUCTION 

Under a uniform gravity, a vertically stratified system of interstellar medium (ISM) composed 
of gas, magnetic field and cosmic-rays is unstable against perturbations (Parker 1966, 1967). This 
nature has been called the Parker instability in the astronomical literature. The instability is 
broken down into two distinct modes; undular and interchange (Hughes & Cattaneo 1987). The 
former, whose wave vectors lie in the plane defined by the direction of gravitational field (the 
vertical direction) and the direction of unperturbed magnetic field (the azimuthal direction), 
undulates the azimuthal field lines and induces the gas to slide down along the field lines into 
“magnetic valleys”. The critical adiabatic index for the undular mode without taking into account 
cosmic-rays is 7cr,u = (1 + «)^/(l + 1.5a), where a is the ratio of magnetic to gas pressure. The 
undular mode has a preferred wavelength along the azimuthal direction. The latter interchange 
mode, whose wave vectors lie in the plane perpendicular to the direction of unperturbed magnetic 
field, leaves the direction of the magnetic field unchanged but alters its strength. The critical 
adiabatic index of the interchange mode is 7cr,i = 1 — a. This mode has a maximum growth rate 
at infinite wavenumber, indicating it is a Rayleigh-Taylor type. If the wave vectors are allowed 
to have all three components, then there exists the mixed mode of the undular and interchange 
modes (Matsumoto et al. 1993). The critical adiabatic index of the mixed mode is 7cr,m = 1 -|- a. 

If three-dimensional (3D) perturbations are applied to a system with the effective adiabatic 
index 7cr,i < 7 < 7cr,u, which is typical of the ISM, the system is unstable to the undular mode but 
stable to the interchange one. The mixed mode, however, has maximum growth rate at vanishing 
wavelength along the radial direction (Parker 1967), showing the characteristic of the interchange 
mode. This implies that the interchange mode as well as the undular mode are involved in the 3D 
development of the Parker instability. It is our understanding that the “magnetic arches”, which 
are formed by the undular mode, are less dense than their surroundings, and are now susceptible 
to the interchange mode. 

It has been considered that the Parker instability is one of the plausible mechanisms for the 
formation of giant molecular clouds (GMCs) (Mouschovias et al. 1974; Blitz & Shu 1980). As the 
formation mechanism of GMCs, the Parker instability, however, has two problems. One is that, as 
we have mentioned in the last paragraph, its linear growth rate is maximum at infinite wavenumber 
along the radial direction of our Galaxy. This means that chaotic structures rather than ordered 
large-scale ones are expected to form. The other is that, according to two-dimensional numerical 
simulations of the Parker instability (Basu et al. 1997), the density enhancement is at most a 
factor of 2 or so, which is too small. 

In this paper, we report the 3D simulations of the Parker instability under a uniform 
gravity, which are intended to study the mixed mode. Our work is a 3D extension of that of 
the two-dimensional undular mode by Basu et al. (1996, 1997). The simulations enable us to 
answer the following questions on the 3D development of the Parker instability; i) What are the 
structures formed in the nonlinear stage? ii) What is the density enhancement factor? Based on 
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the answers, we address the role of the Parker instability in the formation of GMCs. In addition, 
the simulations show that there is a final equilibrium state which is different from that obtained 
in the two-dimensional study (Mouschovias 1974; Basu et al. 1996, 1997). Other numerical works 
on the Parker instability include two-dimensional simulations of the undular mode by Matsumoto 
et al. (1988, 1990) and 3D simulations of the mixed mode by Matsumoto &: Shibata (1992), both 
of which employed a point-mass-dominated gravity. The point-mass-dominated gravity can mimic 
the linearly increasing part of the Galactic gravity at solar neighborhood when the e parameter 
(the ratio of gravitational energy to thermal plus magnetic energy) becomes as large as about 
1000. Their simulations under the gravity correctly demonstrated the nonlinear development of 
the undular and mixed modes, although their numerical resolution was lower than that of Basu 
et al. and ours. However, application of their results to the Galactic disk was hampered, because 
their parameter (e = 6) is too small to model the Galactic gravity. 

In §2, the initial setup and the numerical code are briefly described. The results of the 
simulations and discussions are presented in §3. 


2. NUMERICAL SETUP AND CODE 


To describe the local behavior of the Parker instability, we introduce Cartesian coordinates 
(x, y, z) whose directions are radial, azimuthal and vertical, respectively. Under a uniform external 
gravity, —gz, the initial magnetohydrostatic equilibrium state of a system composed of isothermal 
gas with sound speed, a, and uni-directional magnetic field, Bq = yBo{z), is given by 


Pojz) 

Po(0) 



exp(-z/H), 


( 1 ) 


where the vertical scale height of the gas is dehned as H = (1 -|- a)a^ fg. As Parker (1966) did, we 
assume a{= Bq/ISttpoo^]) to be uniform initially. We study the case with the initial a = 1 and the 
adiabatic index 7 = 1 . 


The computational box has 0 < x < 12ff, 0 < y < 12H and 0 < 2; < 12H, where 12H 
is the azimuthal wavelength of maximum linear growth (Parker 1966). The periodic boundary 
condition is enforced at x = 0, x = 12H, y = 0 and y = 1277, while the reflection boundary 
condition at z = 0 and 2; = 1277. To initiate the instability in the initial equilibrium state, we add 
random velocity perturbations. The standard deviation of each velocity component is set to be 
10“^a. Since the 3D mixed mode has maximum linear growth at infinite wavenumber along the 
radial direction, high resolution simulations are required in order to resolve resulting small scale 
structures. We have performed simulations using 64^, 128^ and 256^ cells, although we report 
mostly the highest resolution one with 256^ cells. Physical quantities of length, speed, density and 
magnetic field are given in units of the density scale height, 77, the isothermal sound speed, a, the 
initial density at 2: = 0, /Oo(0)) and the initial field strength at 2; = 0, 77o(0), respectively. For the 
typical values of a=6.4 km s“^ and 77= 160 pc (Falgarone & Lequeux 1973), the resulting time 



unit, H/a, becomes 2.5 x lO'^ years. The simulations have been made with a 3D isothermal MHD 
code based on the explicit, finite-difference TVD (Total Variation Diminishing) scheme (Kim et 
al. 1998). 


3. RESULTS AND DISCUSSIONS 

Figure 1 shows the time evolution of the logarithmic values of the root-mean-square (rms) 
velocity components in the simulations with 128^ and 256^ cells. Two solid lines are the maximum 
linear growth rates of the mixed mode, Dinax,m = 0.41, and the undular mode, Dinax,u = 0.34 
(Parker 1966, 1967). The evolution patterns in the two panels agree well, indicating the simulations 
produce converged results at least for the global quantities such as The evolution is 

divided into three phases; linear (0 < t ^ 35), nonlinear (35 t ^ 50) and relaxed {t ^ 50). 

During the initial transient stage of the linear phase {0 < t ^ 25), random velocity 
perturbations adjust themselves by generating a one-dimensional magnetosonic wave, which 
moves vertically. The peaks in around t ~ 7 and 21 occur when the wave reaches the 

upper boundary where background density is lowest. After the transient stage, the instability 
grows exponentially (25 ^ t ^ 35). Kim et al. (1998) showed that in a two-dimensional test 
simulation the Parker instability grows at the maximum linear growth rate of the undular mode, 
f^max,u = 0.34. In Figure 1, however, the Parker instability in 3D simulations grows at a rate that 
is smaller than the maximum linear growth rate of the mixed mode, although slightly larger than 
that of the undular mode. This is because the maximum growth of the mixed mode occurs at 
infinite wavenumber along the radial direction, while numerical simulations with finite resolution 
can not resolve such small structures. 

The growth of the instability is slowed down after t ^ 35 and saturated in the nonlinear 
phase due to magnetic tension. The 3D images in Figure 2 depict density structure, velocity 
vectors and magnetic field lines at a nonlinear epoch, t = 40, (Figure 2a) as well as at a relaxed 
epoch, t = 60 (Figure 2b). The two-dimensional images in Figure 3 depict density structure and 
velocity vectors in a sliced yz plane, xz plane, and xy plane at t = 40, when the instability is 
fully developed. Figure 2a shows magnetic arches and valleys created along the magnetic field 
lines. Their azimuthal wavelength is about 12, that of the undular mode. The structures in the 
sliced yz plane of Figure 3 have similarities with those in two-dimensional simulations of the 
undular mode (Basu et al. 1996; Kim et al. 1998). These are the manifestations of the undular 
mode in the 3D evolution of the Parker instability. On the other hand, the iso-density surface in 
Figure 2a is corrugated along the x-direction. Also, in the sliced xz and xy planes, the density 
and velocity changes rapidly along the direction. These are the manifestations of the interchange 
mode. The structure has maximum power at the x-wavenumbers corresponding to 1/16 to 1/32 
of the box size, or 8 to 16 cells. This is the size of the smallest structures the code can resolve 
cleanly (Ryu et al. 1995; Kim et al. 1998). Structures smaller than that tend to be dissipated 
by numerical diffusion. Simulations with higher resolutions would have produced structures with 
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larger x-wavenumbers. Comparison between the results of the high resolution simulation using 
256^ cells and that of the medium resolution simulation using 128^ cells has confirmed the above 
statement. Hence, we can state that due to the manifestation of the interchange mode, the flow 
motion in the 3D simulations should be characterized as chaotic. This is different from what was 
seen in two-dimensional works (Mouschovias 1974; Basu et al. 1996, 1997). 

The chaotic flow motion in the nonlinear phase augments magnetic reconnection, especially 
in the magnetic valley regions (see Jones et al. 1997 for further discussion on reconnection in 
simulations of this kind). As it accumulates, gas in the valleys is compressed and magnetic 
field lines are pushed down by the weight of the compressed gas, generating strongly curved 
field lines. Subsequently, the curved field lines get reconnected, facilitated by the chaotic flow 
motion. The reconnection allows the gas to drop down off the field lines, and at the same time 
the reconnected field lines to float upwards. In this way, gas aggregates along the midplane while 
magnetic field is evacuated from it. The chaotic flow motion induces reconnection in other parts 
of the computational box, but there it contributes less in the redistribution of gas and magnetic 
field. During the relaxed phase, the dense gas with higher pressure than its surroundings spreads 
along the magnetic field lines, and the redistributed gas and magnetic field flatten themselves 
horizontally. A snapshot of the relaxed phase in Figure 2b shows that by t = 60 the corrugation 
has been smeared out, the velocity field has been almost aligned with the magnetic field lines, and 
the magnetic field lines have become more or less straight. However, the iso-density surface, which 
coincides initially with the z = 3 plane, has moved downwards to z ~ 2. This is an indication that 
the relaxed system is significantly different from the initial one. 


To show quantitatively how much gas is segregated from magnetic field, we plot a, which is 
defined as 

<-i 2 H ^ 12 H z;t) 


rizn rl 

a{z]t) = 

Jo Jo 


rl2H r^2H 


-dxdy / 


dxdy, 


( 2 ) 


8Tra^p{x, y, z] t) 

with respect to z at several epochs in Figure 4. Initially a is equal to 1. As time goes by, a 
decreases near the midplane, because gas aggregates there while magnetic field is evacuated as 
stated in the above paragraph. By t = 60, d decreases by a factor of 10 in the midplane. It 
increases more than by a factor of 10 around the upper boundary. At t = 60, in the range 
0 < z ^ 4 magnetic field is more or less uniform and the external gravity is now supported almost 
totally by gas pressure. But for z > 4, ~ 90% of the external gravity is supported by magnetic 
pressure and the rest by gas pressure. 


So, if the Parker instability is allowed to develop into the relaxed phase without being 
interrupted by other events, such as explosive energy injections by OB associations in the Galactic 
plane, the initial unstable Parker system should evolve into a new stable equilibrium state where 
the external gravity is mostly supported by the gas pressure alone. In other words, one role that 
the Parker instability can play in the Galactic plane is to reduce the scale height of interstellar 
clouds while increasing that of interstellar magnetic field. 

Energetics in our simulations has been checked through energy plots (not shown in this paper) 
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which are similar to those in Matsumoto et al. (1988, 1990). The total energy of the system (see 
Mouschovias 1974 for definition) stays constant during the linear phase, but decreases significantly 
during the nonlinear phase. This is because the total energy is no longer a conserved quantity in 
an isothermal flow, when shock waves are generated in the system (Matsumoto et al. 1990). 

How long does it take the Parker instability under the uniform gravity to grow? In estimating 
the elapsed time, the initial transient stage should not be taken into account, since its duration 
depends on initial perturbations and can be reduced if we start simulations with an initial setup 
prescribed by the linear analysis. The density at magnetic valleys starts to increase significantly 
only after the initial transient stage. Hence, the duration of the linear phase except the initial 
transient stage may be regarded as the growth time of the instability in our simulations. It 
amounts to ~ 2.5 x 10® years. What is the density enhancement resulting from the 3D Parker 
instability? At the end of the linear phase, t ~ 35, the density enhancement factor reaches the 
maximum value. Later, it slowly decreases. On the scale of GMCs, ~ 60 pc, which corresponds to 
8 cells in the high resolution simulation, the maximum density enhancement along the midplane is 
^ 2 at t = 35. 

The Parker instability has been considered as one of the plausible mechanisms for the 
formation of GMCs as described in Introduction. Based on our 3D simulations of the Parker 
instability under the uniform gravity, we argue the following: The time scale, ~ 2.5 x 10® years, 
required for the development of the instability is somewhat larger than the lifetime (~ 3 x 10^ 
years) of GMCs (Blitz & Shu 1980). The time scale can, however, be reduced, if we replace 
the uniform gravity by a realistic gravity (Kim & Hong 1998). Therefore, as far as time scale 
is concerned, the Parker instability can be regarded a viable formation mechanism of GMCs. 
However, the structures formed in the 3D simulations are sheet-like or filamentary, and the 
maximum density enhancement factor in the scale of GMCs is only ~ 2. Hence, we conclude that 
it is dijficult to regard the Parker instability alone as the formation mechanism of GMCs. 
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FIGURES 


Fig. 1.— Time evolution of the root-mean-square values of the radial velocity, < the 

azimuthal velocity, < Vy and the vertical velocity, < vl Natural log is used along the 

ordinate. Left panel is from the medium resolution simulation with 128^ cells, and right panel is 
from the high resolution simulation with 256^ cells. Two solid lines represent the maximum linear 
growth rates, 0.34 and 0.41, of the undular mode and the mixed mode, respectively. 

Fig. 2.— Perspective views of density structure, velocity vectors and magnetic field lines at two 
time epochs, (a) t = 40, and (b) t = 60, from the high resolution simulation with 256^ cells. Box is 
oriented in such a way that the radial (x—), the azimuthal (y—), and the vertical (z—) directions 
are from left to right, from near to far, and from bottom to top, respectively. Translucent images 
of density in the three planar slices, xy plane at z = 3, yz plane at x = 6, and zx plane at y = 6, 
are overlapped. Colors are mapped from red to blue as density decreases. Gray surface with an 
iso-density of po(z = 3) is included. Velocity field vectors in the xy plane at z = 3 are represented 
by red arrows. Magnetic field lines, whose starting points lie along the line of z = 6 and y = 0, are 
represented by blue ribbons 

Fig. 3.— Images of density structure and velocity vectors on the three planes, x = 6 (left), y = 6 
(center), and z = 3 (right) at t = 40 from the high resolution simulation with 256^ cells. The three 
planes are same as the sliced planes in Figure 2. Colors are mapped from red to blue as density 
decreases. The unit of the velocity vectors is shown at the right lower corner of the image. 

Fig. 4.— The ratio of magnetic to gas pressure a averaged over xy plane (see eq. i) at several 
time epochs from the high resolution simulation with 256^ cells. Common log is used along the 
ordinate. 
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This figure "fig2.gif" is available in "gif" format from: 


http://arXiv.org/ps/astro-ph/9808244vl 


























